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Abstract 


A 1D dynamic solid oxide fuel cell (SOFC) model has been developed for real time applications. The model accounts for all transport and 
polarization phenomena by developing a system of governing differential equations over 1D control volumes. The 1D model is an improvement 
over existing OD real time models in that it can more accurately predict the temperature and pressure variations along the cell while maintaining 
real time capabilities with regards to computational time. Several simplifications are required to maintain real time capabilities while improving 
the fidelity of the model. 

It was found that a 1D model with 21 nodes performs each time dependent computation in 3.8 ms. Results show that activation overpotentials 
account for most of the total cell overpotential, and that temperature variations across this particular cell exceed 100 K. Depending on the gas 
channel configuration, the pumping power required to supply air and fuel to the cell can be in the same order of magnitude as the power produced 
by the fuel cell. It is shown that the use of fewer channels with larger cross-sectional areas, reduces the pressure drop due to wall friction, and 


hence reduces the required pumping power. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Solid oxide fuel cells (SOFC) are considered prime can- 
didates for stationary power generation in the intermediate to 
long-term future. In addition to its clean and efficient operation, 
its high temperature of operation (800-1100 K) allows for use 
of a wide range of fuels (including CO), as well as the possi- 
bility of heat recovery and co-generation. Because of its high 
tolerance to impurities in the fuel, e.g. CO and S, as compared 
to polymer electrolyte membrane (PEM) fuel cells—SOFCs can 
also be used in hybrid power generating systems. It can be cou- 
pled with a gas turbine or a biomass gasifier, and as such can be 
integrated with other renewable technologies. 

Fuel cell modeling plays a critical role in the development 
of fuel cell technology. Fuel cell models can be broadly catego- 
rized as transport models or system models. Transport models 
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provide a more detailed analysis of the transport phenomena 
occurring within the cell and can be used to optimize cell design 
and operating conditions. System models simulate the behavior 
of the fuel cell within the overall power generating system. 

System models are not as comprehensive as transport models 
because they do not require the details of the cell’s internal phe- 
nomena, only the input/output characteristics. Transport models 
can range from 1D to 3D and may by steady state or transient 
[1-8]. System models are usually OD (lumped) or 1D to mini- 
mize the memory and computational time requirements, and are 
almost always transient [8—11]. In a system model, it is usually 
required that the computations be performed hundreds of times, 
therefore it is required to keep the memory requirements and 
computational times to a minimum. Further discussion of trans- 
port and system modeling is given in Bove and Ubertini [12] 
and Liese et al. [13]. 

Real time modeling is a special category of system modeling 
where the fuel cell model operates in real time, while fully inter- 
acting with other equipment within the system. This requires 
that the time dependent computations be performed within a 
specified period of time—typically in the order of milliseconds. 
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Nomenclature 


Arabic symbols 

area (m?) 

specific heat capacity (Jkg~! K7!) 
diffusivity (m? s7!) 

Activation Energy (J mol7!) 

friction factor 

Faraday constant (96485 J mol~!) 

heat transfer co-efficient (W m~? K7!) 
cell current (A) 

thermal conductivity (W m`! K7!) 
exchange current density term (S m7?) 
length (m) 

mass flow rate (kg s7!) 

number 

pressure (Pa) 

heat generation (W) 

universal gas constant (8.3143 J mol`! K7!) 
Reynolds number 

time (s) 

temperature (K) 

velocity (m s7!) 

potential (V) 

width (m) 

mole fraction, displacement in the direction of the 
gas channel (m) 


DOTI DASEN Np > 


Bo) 
A 


sexes 


Greek symbols 

À stoichiometric rate 

u dynamic viscosity (Pa s) 
p density (kg m7?) 

T thickness (m) 

w mass fraction 


Subscripts and superscripts 


act activation 

an anode 

atm atmospheric 
ca cathode 

ch channel 

el electrode 

gen generated 

i species 

int interconnect 
ohm ohmic 

pen positive-electrolyte-negative 
0 reference state 


Most system models are designed to perform computations in 
the order of seconds and minutes. In fact most system mod- 
els are designed to minimize memory requirements rather than 
maximize speed. The precise time frame which defines “real 
time” depends on the other equipment in a fuel cell system. In 
this work, it depends on the compressor valve timing, which is 
approximately 5 ms. 


The hardware-in-the-loop (HIL) concept makes use of real 
and virtual components of the fuel cell system. In the example 
of a turbine/SOFC hybrid, the turbine power plant comprises 
the real components. The real time fuel cell model comprises 
the virtual components. Thus the model operates in lieu of a 
real fuel cell stack, and interacts in real time with the other 
components. The other components operate as if there were, in 
fact, an actual fuel cell stack present. Thus, for HIL purposes, it 
is critical that the model operates in real time. The fact the real 
time model must fully interact with other equipment means that 
cell voltage and current are not the only input/output parameters 
of interest. Also important are input and output temperatures, 
pressures, concentrations and flow rates. The output from the 
model is inputted into real components, which then perform 
as though there were in fact a fuel cell present outputting at a 
given temperature, pressure, concentration and flow rate. So in 
addition to current/voltage, the real time model is also interested 
in various transport properties. 

It is believed that hybrid fuel cell systems will play a promi- 
nent role in the intermediate future as fuel cell technology 
slowly makes its way into the market, and real time modeling 
will become increasingly important. Most of the system models 
developed are designed to deliver “fast” computations, but not 
fast enough for real time simulation. Most system models per- 
form the necessary computations in the order of seconds. Real 
time models, therefore, must be simplified to allow for compu- 
tational times in the order of milliseconds. For this reason, most 
real time models are OD [14,15]. The most popular hybrid system 
investigated is the SOFC/gas turbine system [16,17], however 
biomass/SOFC systems have also received some attention [18]. 
In the literature, real time models for hybrid power generating 
systems have focused on control systems [19] and have covered 
a wide range of applications including aerospace applications 
[20]. 

This work presents a 1D real time model, which is designed 
to improve the fidelity of real time modeling, while maintaining 
the necessary computational times. As mentioned earlier, the 
fidelity of the model must be measured in terms of its ability to 
predict the electrochemical as well as the transport performance 
of the fuel cell. Results are validated by comparing the prediction 
of the real time model to those of a fully resolved model using 
a finite element mesh with iterative solvers. The model used to 
perform the validation is much more comprehensive and does 
not make many of the limiting assumptions of the real time 
model discussed in this paper. Thus, we are validating how well 
the real time model compares with a more comprehensive model 
of the same problem. The model is presented here for a planar 
cell, however, it can be applied to any type of SOFC design. 


2. Model development 
2.1. Assumptions and simplifications 

Compared to other system models, real time models must be 
simplified in order to perform the computations within the nec- 


essary time frame. However, compared to existing OD real time 
models, the present 1D model is an improvement with regard 
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to accuracy, while maintaining real time capability. This section 
discusses some of the simplifications employed in the present 
model. 

The gas channel and interconnect are integrated at both the 
fuel and oxidant side. The treatment is still 1D, but the chan- 
nel and interconnect states are not separated. For example, the 
results presented by Aguiar et al. [7] show a difference between 
the interconnect temperature and the fluid temperature in the 
channels. This work treats the channel/interconnect as one inte- 
grated domain. 

The anode, cathode and electrolyte are similarly integrated 
into one region. This region is called the PEN. For the purpose 
of calculating pressure drops in the porous electrode regions, the 
thicknesses of the anode and cathode are used. However, there 
is no distinction made between the electrolyte temperature and 
the electrode temperature. 

To reduce the number of computations, a uniform distribu- 
tion of current density is assumed. The average current density 
is assumed to exist at all points. This way there is no need to 
iteratively calculate the cell current or voltage. This is in con- 
trast to numerous models which assume equipotential conditions 
with a spatial variation in current density. These models, hence 
will be more accurate, but it will be seen in Section 3 that the 
inaccuracy of the present model is not severely compromised by 
this assumption. On the other hand its computational speed is 
greatly enhanced. Multi-dimensional models, and even the 1D 
model by Aguiar et al. [7] require an iterative solver to con- 
verge upon the correct solution to the electrochemical problem. 
This is not amenable to real time simulation, since the iterative 
procedure is the primary time consumer in the solution process. 

For the purpose of the present model, it is not necessary to 
know the overpotentials at each node, only the average overpo- 
tentials across the cell. The reversible potential as well as the 
activation and ohmic overpotentials depend on temperature and 
partial pressure. The temperature and pressure will be computed 
at each node, however, to determine the necessary potentials 
and overpotentials, the average PEN temperature and partial 
pressures will be used. 

Following the previous two simplifications, the heat genera- 
tion due to cell irreversibilities is also assumed to be uniformly 
distributed and equal at each nodal point. Also uniform heat and 
mass transfers between the PEN and the channel is assumed at 
each node. 

Due to the relative time scales, the electrochemical problem 
is assumed to occur instantaneously (or in steady state). It is 
assumed that the voltage immediately responds to changes in 
current. However, the thermal problem takes a longer period of 
time to reach steady state conditions, and the electrochemical 
problem depends on temperature. So the potentials are consid- 
ered to depend on time only as the temperature depends on time. 


2.2. Electrochemical model 


The following half-cell reactions take place at the anode and 
cathode, respectively: 


2H? +20% > 4e~ + 2H20 (1) 


O2 +4e7 > 2077 (2) 


In this model, the cell current is specified and the cell over- 
potentials are then determined. The reversible cell potential can 
be determined from Eqs. (3) and (4) [14]: 


RT, ( Pa, Po; 

E = Eo + —In -a (3) 
2F  \ Pro Patin 

Eo = 1.2877 — (2.904 x 10~4)T (4) 


The average channel partial pressures are used in Eq. (3). The 
PEN temperature is used since that is where the electrochemical 
reactions occur. To determine the concentration overpotential, 
the average partial pressures in the PEN can be substituted in 
Eq. (3). The concentration overpotential is taken as the difference 
in values obtained by using the channel pressures and the PEN 
pressures in Eq. (3). 

The average PEN partial pressure of species, i, can be deter- 
mined from its average partial pressure in the channel using 
Fick’s law of diffusion, assuming that convection is negligible 
within the porous electrodes. Previous works consider the gas 
concentration as a function of partial pressure only [7,8]. How- 
ever, based on the ideal law, gas concentration is a function of 
partial pressure and temperature: 

oO 

Mi T 5) 
t M; Di ( 


pee 
Rig RTs, 


The ohmic overpotential is a function of current density and 
temperature. Eq. (6) includes the PEN resistance as well as the 
interconnect resistance [7,21]: 


T 


10.3 x 10° 
Nohm = i tox x 107!) exp (=) 


+ (0.23174) es(-001157) (6) 


The expressions for the activation overpotentials are taken 
from Liese et al. [14]. It is assumed that there is no partial pres- 
sure dependence at the cathode. The exchange current density 
is given by Eq. (9). In Eqs. (7) and (8), the PEN partial pres- 
sures are used. The PEN subscript is not shown to avoid over 
complicating the equations: 


an | i/2io 
Nact = sinh ref ref 
F (Pa, / PR, )(Pmo/ PR o) 


1 P pref 
in { 7! a (7) 
2 Pm o/ Piho 
© RTĪ. a 
naa = -F sinh (=-)| (8) 
. RT, Eel (9) 
= — Kaexp | -— 
n= aF oe RE 
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Fig. 1. Energy and mass balances over 1D control volumes. 


2.3. Transport model 


The transport equations are developed based on mass and 
energy balances over the 1D nodal control volumes shown in 
Fig. 1. Because there is mass transfer between the channels and 
the PEN, there is a non-zero source term in the continuity equa- 
tion (10). This source term is needed to account for mass transfer 
in the orthogonal directions: 


aou) _ 9 Ymi 
Ox ot = AchLcel 


(10) 


The sign convention used in this work is such that a vector 
quantity is considered positive when the direction is into the 
domain, and negative out of the domain. The rates of diffusion 
of each individual species in the PEN are as follows: 


= 11 
mH, o ÍT (1) 
ò IMu, 
= — = 12 
mH, 2F ( ) 
o IMo, 
= 2 13 
O2 4F ( ) 


The diffusion of hydrogen and oxygen are considered nega- 
tive because they flow out of the channel domains, while water 
vapor flows into the anode channel domain: 


The momentum equation (14) is based on the Navier-Stokes 
equations applied over 1D. Wall friction is based on the velocity 
gradients in the directions perpendicular to the solid surfaces, 
thus the viscous terms must be modified to properly account 
for wall friction in the channels. It can be shown that for a 
rectangular channel, the last term in Eq. (14) accounts for the 
pressure drop due to Darcy friction, assuming laminar flow. Eq. 
(15) gives the dimensionless product of the friction factor and 
Reynold’s number based on analytical correlations given in [22]. 
The Reynold’s number is based on the square root of the channel 
cross sectional area. Note that the Aen in Eq. (14) refers to the 
total cross-sectional area of all of the channels, so that equation 
is modified by multiplying by the total number of channels. 

The energy Eqs. (16) and (17) govern the thermal problem. 
All of the heat generation is assumed to occur in the PEN, where 
conduction is dominant and convection is considered negligible. 
The channel/interconnect cross-sectional area ratio occurring in 
the first term on the RHS of Eq. (16) arises due to the inte- 
gration of the channel and interconnect domains. It can also be 
seen from these equations that heat conduction is considered 
negligible in the fluid phase, due to the fact that the thermal 
conductivity is much higher in the solid regions than in the fluid 
phases. The same is true of the heat storage due to higher vol- 
umetric heat capacities (pcp) of the solid regions. Also the sign 
convention of positive heat and mass transport into the domains 
are maintained: 


2 o (pcp). OTint = O(pucTint) Ach +k 8? Fant 
aP A Apu? Ehu Din oe A, ta 
ox ot ox Ach Leell ò 
1/2 1/2 X mci Tint + hAcel(Tpen — Tint) 16 
_ | teh í a) ast) : f Re HNeb (14) * AintL cell ii 
Weh Tch Ce Aas 
12 
S Re Kawa, = (15) 


[1 — (192/75)(Ten/ Wen) tanh((7/2)(weh/Teh))] CL + (Toh/ Weh))V Teh/Weh 
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OT pen 3” Tpen 
(Pcp)pen ~a; = Kena 9 
o 
+ Qegen = > mici Tin — hAcet(Tpen — Tint) 


Apen Leet 
a7) 


The heat generation in the cell can be determined from first 
and second law thermodynamic analyses of the fuel cell. The 
heat is due to irreversibilities due to the total cell overpotential 
and entropy changes. The entropy change for the fuel cell reac- 
tion can be correlated from thermodynamic tables and expressed 
as a function of temperature. This expression is given in Eq. (19): 


=I V, V, Eee (18) 
Q gen = rev — Vcell — AF 
AS -ly-l 
=" = —(23.328 + 0.0042T) (J mol`! K7!) (19) 
n 


The Stefan-Maxwell equation (20) is used to describe the 
transport of individual species in the gas channels. The species 
considered are H2, CO2 and H20 at the anode; and O7 and N2 at 
the cathode. In this paper, it is assumed that the fuel is externally 
reformed resulting in an inlet fuel stream consisting of H2, CO2 
(inert) and H20. No direct internal reforming is considered in 
this paper, as that will be addressed in a later work. 

The last term in Eq. (20) is the source/sink term which rep- 
resents mass diffusion between the channels and the PEN. No 
electrochemical reactions take place in the channels, but it is 
assumed that all of the net mass of a species diffusing from the 
channel to the PEN eventually reacts: 


dwi puw) Cowi) Mi 


Py T ax ' 8x2 AchLeell 


(20) 


Finally the gas mixture density is given by the ideal gas 
relation: 


PM P oN T! 
m z 21 
P= RTO RT (£ e) 21) 


2.4. Boundary conditions 


The boundary conditions at node 1 are specified from the 
inlet gas conditions. The anode and cathode gas temperatures, 
mass fractions, pressures and flow rates are given (Dirichlet). 
The PEN temperature at node | is not specified, but assuming 
negligible heat losses at the ends of the cell, insulation conditions 
exist (Neumann). 

The same is true of the PEN at the final node. For the chan- 
nel/interconnect, it is assumed that heat and mass exit the cell 
solely via convection, i.e. heat conduction and mass diffusion 
are zero at the outlet (Neumann conditions). 

Tables 1-3 list the numerical values used for base calculations 
in the present model. These numerical values are for the most 
part the same as those used by Aguiar et al. [7]. 


Table 1 
Cell geometries [7] 


Number of channels per cell 50 

Cell length 0.4m 
Cell width 0.1m 
Channel thickness 0.5mm 
Channel width 1mm 
Interconnect thickness 2.5 mm 
Anode thickness 500 um 
Cathode thickness 50 um 
Electrolyte thickness 20 um 


Table 2 
Cell constants [7] 


Tin K) 

Pin (atm) 

Ka (Sm~) 

Ee (kJ mol!) 
h 

h (Wm? KT!) 
D (m? s7!) 


Anode 


1023 
1 
6.54 x 10!! 
140 
1.33 
302 
3.66 x 1075 


Cathode 


1023 
1 
2.35 x 10!! 
137 
8.93 
302 
1.37 x 1075 


Interconnect PEN 


pcp (MJ m~? K7!) 4 2.95 
k (Wm! K7!) 25 2 


2.5. Solution methodology 


It is noticed from Eqs. (10) and (14) that the transient term in 
the continuity and momentum equations are placed on the RHS 
rather than the LHS. This was done deliberately to allow for 
greater stability in the solution process. The momentum equa- 
tion (14), for example, has two dominant terms—the pressure 
gradient and the Darcy friction. Placing the transient term on 
the LHS requires very precise initial guesses, which place severe 
limitations on the model. There would also be numerical stability 
problems when rapid changes in loading conditions are applied. 
Physically, it is best to understand the momentum equation as 
a formula for calculating the pressure drop given the velocity 
profile. 

Recall that in the present model, iterative solution procedures 
are avoided. Iterative procedures require time to converge, and 
are thus not amenable to real time modeling. An explicit solu- 
tion scheme is preferred, which allows the model to calculate 
the next time parameters given the present state variables and 
derivatives. It is assumed that the necessary history of the cell 


Table 3 
Gas properties [23] 
Inlet mole Specific heat capacity, Dynamic 
fraction, x cp I kg! K7!) viscosity, u (Pas) 
H2 0.571 13960 + 0.950T 5.24+0.0156T 
H20 0.286 1639.2 +0.641T —9.56 + 0.3603T 
CO2 0.143 804+ 0.460T 70.40 + 0.3340T 
On 0.21 876.80.217T 146.38 + 0.3330T 
N2 0.79 935.6 +0.232T 129.22 + 0.2725T 
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Fig. 2. Solution algorithm. 


can be captured by its present state variables and time deriva- 
tives. The solution algorithm is shown in Fig. 2. The model 
is solved in Matlab/Simulink using the 4th order Runge-Kutta 
method to solve the time dependent problem. The fully resolved 
model used to validation of the real time model is generated for 
a finite element mesh, and solved iteratively using FEMLAB® 
3.1. All computations are conducted on a 1GB, 3GHz Pentium 
4 processor. 


3. Results 
3.1. Computational time 


Fig. 3 shows the variation of the steady state channel outlet 
temperatures as the number of nodal points in the model varies. 
As expected, as the number of nodes increases, the solution 
converges to the point of grid independence. It shows that using 
less than 10 nodes results in large inaccuracies. In this range, the 
solution is very grid dependent. However, as the number of nodes 


increases above 20, the solution becomes less grid dependent. In 
fact, there is very little difference in accuracy between a model 
with 21 nodes and one with 51 nodes. A model with 2 or 3 
nodes approaches a OD model. These results clearly illustrate 
that a 1D model greatly improves the accuracy compared with 
a OD model. 

Fig. 3 also shows the computational times required to perform 
each time dependent iteration. This value increases approxi- 
mately to the 1.5 power of the number of nodes. Given that 
5 ms is the desired maximum computational time per iteration, 
the model with 21 nodes provides the optimum balance between 
speed (3.8 ms) and accuracy. Its accuracy is not very different 
from one with 51 nodes, but its computational time is signif- 
icantly lower. Thus a 1D model with 21 nodes improves the 
accuracy of OD models while maintaining real time capabilities, 
and thus 21 is chosen as the optimum number of nodes and will 
be used for the remainder of results presented. Note that the 
number of nodes is chosen as 21 so that there will be 20 spatial 
divisions. 


3.2. Model validation 


Validation of the time dependent model is done by com- 
paring numerical results with the results of a fully resolved 
model for the same set of operating conditions and cell geome- 
tries. This fully resolved model uses an iterative technique to 
solve for the current density distribution for a specified set of 
equipotential condition. This is done to ensure that the limiting 
assumptions made in the real time model did not significantly 
affect the results. Figs. 4 and 5 show that the inaccuracies intro- 
duced by this assumption only resulted in an error in the order 
of 1%. However, the more comprehensive model requires com- 
putational times between 2 and 3s for a 1D domain. This far 
exceeds the requirements for real time modeling. These figures 
show that the real time model’s accuracy is not significantly 
compromised by the simplifying assumptions. 

Fig. 4 shows the dynamic voltage response to changes in cell 
current. Initially, steady state conditions were allowed for a cell 
current density of 0.5 A cm~? (cell current of 200 A). Then the 
current density is suddenly changed to each of the other values 
shown in the figure. Consistent with model assumptions, when 
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Fig. 3. Solution accuracy and computational times vs. the number of nodes. 
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Fig. 4. Dynamic cell potential response to change in loading conditions (constant 
heat capacity). 


the current density changes, the voltage immediately responds. 
However, because the temperature takes a longer period of time 
to reach steady state, the voltage slowly stabilizes as the tem- 
perature varies with time. Thus, there is the expected overshoot 
and undershoot of the cell voltage followed by a period of sta- 
bilization. 


3.3. Thermal analyses 


The results in Fig. 4 are based on the assumption that the 
specific heat capacity of the gases do not change significantly 
over the length of the cell—an assumption made by Aguiar et 
al. [7]. So Fig. 4 was produced using the heat capacity values 
based on the inlet temperatures just to compare with their model. 
In reality, heat capacity is a function of temperature, and the 
temperature variation across the length of the cell can be in the 
order of 100K. Table 3 shows the heat capacity correlations 
used in this work. For the case where “constant heat capacity” 
is used, the specific heat capacity was mass averaged at both 
the anode and cathode inlets. In the case where “variable heat 
capacity” was used, the heat capacity is mass averaged at every 
node, taking into account composition and temperature effects. 
Fig. 5 shows the transient voltage response to the same changes 
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Fig. 5. Dynamic cell potential response to load changes (heat capacity as a 
function of temperature and composition). 
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Fig. 6. Dynamic anode outlet temperature response to load changes. 
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Fig. 7. Dynamic cathode outlet temperature response to load changes. 


in current density as in Fig. 4, this time with gas species heat 
capacity as functions of temperature. The overall mixture heat 
capacity is a function of temperature and gas composition. It can 
be seen that the voltages are slightly higher than those in Fig. 4. 
The reason for this is that the temperature distribution is altered 
by the assumption of variable heat capacity, hence the potentials 
and overpotentials are also altered. 
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g. 8. Dynamic PEN (maximum) temperature response to load changes. 
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Figs. 6—8 show the dynamic temperature response to the iden- 
tical changes in current density, for both cases of constant and 
variable heat capacity. They all show that the outlet tempera- 
tures of the channels as well as the maximum PEN temperature 
increase when variable heat capacities are taken into account. 
From the heat capacity correlations in Table 3, all of the gas 
heat capacities increase with temperature. One would expect 
that if the heat capacity increases, then the outlet temperatures 
should decrease. The reason the temperatures actually increase is 
because the mixture heat capacity of the anode decreases consid- 
erably along the channel. Since hydrogen is being consumed, and 
the hydrogen heat capacity is much higher than the other gases, 
the decreasing hydrogen mass fraction results in a decreasing 
mixture heat capacity along the anode channel. Fig. 9 shows the 
spatial variation in mixture heat capacity along the gas chan- 
nels. The cathode channel heat capacity only slightly increases 
along the channel, whereas the anode heat capacity decreases 
significantly as the hydrogen composition diminishes. Since the 
temperatures at both channels and the PEN are all interrelated, 
the net effect is that all of the temperature variables increase 
when the heat capacity is considered a function of temperature. 
One would expect that this consideration would improve the 
accuracy of the model. 

It is noticeable, comparing Figs. 6 and 7, that the anode chan- 
nel outlet has a higher temperature than the cathode. It is true, 
for a SOFC that the reaction heat is more concentrated on the 
anode side, however the present model, which integrates the 
entire PEN region, does not “know” this. The reason the anode 
has a higher temperature than the cathode is because of the lower 
flow rates used there. Table 4 shows that the volumetric flow rates 
at the cathode are an order of magnitude higher than those at the 
anode. This is because the stoichiometric rates are higher at the 
cathode. The anode stoichiometry is chosen for optimum fuel 
utilization, while the cathode stoichiometry is chosen to reduce 
the temperature. 

Another noticeable feature of Fig. 7 is the kink in the temper- 
ature variation at the instant the current density changes. This 
kink exists because the cell current is assumed to change instan- 
taneously, and as a result the rate of mass diffusion between the 
channels and the PEN also changes instantly. When the current 
suddenly increases, the rate at which the cathode channel loses 
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Fig. 9. Variation of steady state mixture heat capacity along the gas channels. 


Table 4 

Parametric model outputs for various input conditions 

Inlet temperature (K) 1023 (base) 973 1073 1023 

Inlet pressure (mV) 1 1 1 3.6 

Cell power (mV) 137.4 120.4 149.2 153.6 

Cell potential (mV) 687.0 602.0 746.0 767.9 

Reversible potential (mV) 918.4 931.8 903.9 950.2 

Activation overpotential 180.8 260.7 119.9 129.5 
(mV) 

Ohmic overpotential (mV) 43.1 62.2 29.8 45.4 

Concentration overpotential Ta 6.9 8.2 T3 
(mV) 

Anode outlet temperature 1144.9 1114.5 1181.0 1134.3 
(K) 

Cathode outlet temperature 1137.3 1105.9 1174.0 1127.3 
(K) 

Maximum PEN temperature 1147.9 1117.5 1183.9 1137.1 
(K) 

Heat generation (W) 109.2 126.2 97.4 99.0 

Anode heat gain (W) 28.3 29.5 27.9 26.9 

Cathode heat gain (W) 81.6 97.5 70.1 72.7 

Anode pressure drop (kPa) 8.6 7.9 9.3 2.4 

Cathode pressure drop (kPa) 92.1 85.7 98.9 25.5 

Anode inlet volumetric flow 203.1 193.1 213.0 56.4 
rate (cm? s7!) 

Cathode inlet volumetric 1850.3 1759.8 1940.7 514.0 
flow rate (cm? s7!) 

Required pumping power 172.1 152.3 194.0 13.2 


(W) 


energy instantly increases, which results in an initial tendency 
for the temperature to drop. However, because more heat is pro- 
duced in the cell when the current increases, the long-term trend 
is for the temperature to increase. The kink represents a conflict 
between the initial loss of energy and the long-term increase in 
heat generation. 


3.4. Parametric analyses 


Table 4 shows the breakdown of steady state cell parameters 
for various inlet conditions. The second column (inlet tempera- 
ture of 1023 K) represents the base case condition. For a current 
density of 0.5 A cm7?, the cell potential is 0.687 V resulting in 
137.4 W of electrical power and 109.2 W of heat generation. The 
overpotentials are dominated by the activation overpotential, 
which constitutes 0.181 V. 

As the inlet temperatures increase, the reversible potential 
decreases, as one would expect from Eq. (4). However, the cell 
potential increases because the activation and ohmic overpoten- 
tials decrease with temperature. The concentration overpotential 
increases slightly with temperature because, according to Eq. 
(5), the reactant concentrations decrease slightly as the temper- 
ature increases. When the inlet pressure increases, the reversible 
potential predictably increases, while the activation overpoten- 
tial decreases, resulting in a net increase in cell voltage. 

The computations for the base conditions show that the heat 
gains by the anode and cathode are 28.3 and 81.6 W, respec- 
tively, representing a 0.6% error in heat balance. The anode 
absorbs less heat than the cathode, yet attains a higher outlet 
temperature. The reason for this, as previously explained, is 
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Table 5 
Effect of channel configuration 


Number of channels Channel thickness (mm) Channel width (mm) fRe [16] Required pumping power (W) 
50 0.5 1 16.457 172.1 
50 0.75 1 14.568 64.5 
50 1 1 14.132 34.8 
25 1 2 16.457 21.5 
10 1 5 25.569 17.0 
5 1 10 36.807 15.8 


the low flow rates existing at the anode. When the inlet tem- 
peratures are increased, there is an increase in heat generation 
primarily because of the high irreversibilities due to entropy 
changes. When the inlet pressure increases from | to 3.6 atm, 
there is a decrease in heat generation and outlet temperatures. 
This is because the better cell performance results in a higher cell 
potential and lower overpotentials. The lower total overpoten- 
tial results in the lower heat generation value, and consequently, 
lower outlet temperatures. 


3.5. Pumping power 


The pressure drop along the channels and the required pump- 
ing power to overcome this pressure drop have not been given 
much attention in previous works. However given the numbers in 
Table 4, this represents a major loss. For the channel geometries 
given, the computed pressure drop across the cathode is a stag- 
gering 92.1 kPa for the base condition. The 8.6 kPa drop across 
the anode is because of the lower flow rate. For the necessary 
volumetric flow rates, the total pumping power required to sup- 
ply the anode and cathode gases to the cell is 172.1 W per cell. 
Considering that the cell power output is 137.4 W, the given cell 
uses more energy than it produces. The other results show that 
the pressure drop increases with temperature, perhaps because 
higher temperatures result in lower gas densities, which result 
in higher velocities, and hence larger wall friction losses. When 
the inlet pressure in raised to 3.6 atm, the pressure drops and 
pumping power decrease significantly. The reasoning is similar. 
A higher pressure results in a higher gas density, which eventu- 
ally reduces the effect of wall friction. Appendix A shows some 
rough “hand” calculations to verify the model results, and to 
show that for the given geometries and operating conditions, the 
pressure drops shown in Table 4 are indeed expected. 

Clearly such a high amount of required pumping power is 
unacceptable. Table 5 shows the effect of configuring the chan- 
nel geometry on the required pumping power. For each case, the 
cell power is the same as the base condition, 137.4 W. The top 
row represents the base condition, where there are 50 channels of 
2 x 1 rectangular sections. For the base condition, as mentioned 
previously, the required pumping power is 172.1 W. Keeping 
the channel width fixed and increasing its thickness results in a 
decrease in pumping requirements. It is primarily because the 
larger channel area allows for lower velocities, hence lower pres- 
sure drops. Now keeping the channel thickness fixed at 1 mm 
and increasing the channel width (hence decreasing the number 
of channels), results in a further decrease in pumping require- 


ments. Although the total cross-sectional area remains the same, 
fewer channels means fewer walls, hence lower pressure losses. 
These results indicate that as far as pumping requirements are 
concerned, having fewer thicker channels is the most optimal 
configuration. 

The results presented in Table 5 do not take into account 
changes in the convective heat transfer coefficient incurred by 
the changes in geometry. It is known that the Nusselt num- 
ber may vary by up to 100% over the geometries mentioned in 
Table 5. However, this does not significantly affect the temper- 
ature distribution. The heat generation is unaffected by changes 
in geometry, thus the only effect will be to minimize the dif- 
ference between the interconnect/channel temperature and the 
PEN temperature. For the geometry shown in the bottom row of 
Table 5 (1 x 10 channel), the steady state temperatures predicted 
by the present model are: anode outlet (1145.8 K), cathode outlet 
(1137.7 K), PEN node 21 (1148.6 K). With a slight modification 
in the script, the effect of channel geometry on convective heat 
transfer coefficient can be taken into account. The respective 
steady state results are: anode outlet (1143.6 K), cathode outlet 
(1139.4 K), PEN node 21 (1144.9 K). It is clear that the effect of 
this phenomenon the temperature distribution is not significant. 
The difference between the three temperatures decreased, but 
not significantly. 


4. Conclusion 


A 1D SOFC model was developed for application in real time 
system simulation. The 1D model is an improvement over exist- 
ing OD real time models. Results showed that a model with 21 
nodes offered the best balance between computational time and 
accuracy. Such a model required 3.8 ms of computational time 
for each time iteration, and offered grid independent solutions. 

The model applied the principles of mass and energy conser- 
vation over 1D control volumes to arrive at a system of governing 
differential equations which were solved in Matlab/Simulink. 
Thus, the model was able to more accurately predict the transport 
phenomena than existing real time models. The model results 
showed that the expected temperature variations for the fuel cell 
in question are in the order of 100 K. Depending on the channel 
geometry, the required pumping power to overcome Darcy wall 
friction was in the same order of magnitude as the power pro- 
duced by the cell. The model predicted that using fewer channels 
with larger cross sectional areas would significantly reduce the 
pressure drop and hence increase the net power produced by the 
cell. 
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Appendix A 


Rough calculations for the pressure drop and pumping 
requirements: 


Channel thickness: Ten = 0.5e — 3 m 


Channel width: Weh = le — 3 m 


f Re = 16.457 


Characteristic length for Reynold’s number : 


d = J/ThWeh = 0.71e — 3m 


At the cathode, for a cell current of 200 A, an inlet O2 mass 
fraction of 0.22, and a stoichiometric rate of 8.93, the inlet flow 
rate at the cathode is 


(200) x (32e — 3) 1 
4 x (96,485) 0.22 
(1.013e5) x (29.16e — 3) 
= (8.3143) x (1023) 
6.73e — 4 
" = (0.347) x (0.5e — 6) x (50) 
u 7 4e — 5Pas, 
Ran (0.347) x (77.5) x (0.71e — 3) 
4e—5 
which verifies that the flow is laminar: 


16.457 
f~ —— = 0.035 
476 


puAch = 8.93 x = 6.73e — 4kgs"!, 


= 0.347kgm >, 


=77.5ms_!, 


= 476 


The wall friction shear stress : 


(0.347) x (17.5) 


Owall © 0.035 = 36.0 Pa 


Friction force per channel length 


= (36.0) x (le — 3 + 0.5e — 3) x 2 = 0.108 N m7! 


Pressure drop per channel length 


0.108 


= = 216, 058 Pam! 
(le — 3) x (0.5e — 3) = 


Over a length of 0.4m, the pressure drop=86,423 Pa 
= 86.42 kPa. 


This is in the same order of magnitude as that predicted in 
Table 4. 


The volumetric flow rate 


= (77.5) x (0.5e — 6) x (50) = 1.94e — 3 m? s7! 


Thus the required pumping power for the cathode 


= (1.94e — 3) x (86423) = 167.8 W 


These calculations show that the pressure drops and pumping 
requirements predicted by the model are accurate for the given 
cell geometries and operating conditions. 
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